Demographic noise as a triggering mechanism for algal blooms 
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Population models, such as those for plankton dynamics, are often based on a mean-field ap- 
proximation of individual behaviors. A weakly stable mean-field configuration, however, can be 
destabilized by demographic noise. In certain cases, such destabilization persists even in the ther- 
modynamic limit. The possibile relevance of such an effect to the onset of algal blooms is discussed. 
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I. INTRODUCTION 

Phytoplankton provide the basis of the food chain in 
the ocean. They are constituted for the large part by 
microscopic algae and their abundance is such that they 
are responsible for roughly one half of the total photo- 
synthesis going on in the planet [1]. 

Although phytoplankton are present in virtually all of 
the world hydrosphere (more precisely, the top 100 m 
layer of the water column, called the euphotic layer, 
where there is enough light for photosynthesis), their dis- 
tribution in space and time is far from uniform [2] . Space 
patchiness is observed at scales that go from that of the 
individual to several hundreds kilometers [3, 4]. The spa- 
tial patterns, revealed e.g. by remote sensing, include 
filaments, fronts and more irregular shapes, suggesting 
that turbulent transport by sea currents may play an 
important role [5-9]. 

Variability in time, on the other hand, is characterized 
by sporadic bloom events in which the plankton density 
in extended regions can grow by up to three orders of 
magnitude in few days [10]. These are seasonal events 
that may or may not recur annually, and are superim- 
posed to the weaker day-night and seasonal cycles. Their 
importance is utmost for several reasons. Depending on 
the species (e.g. diatoms), due to their abundance, bloom 
events can have consequences even at the level of the 
global biogeochemical cycle [11]. Bloom events involving 
other species (e.g. dinoflagellates) can have harmful ef- 
fects on water quality and fishery [12, 13]. Being able 
to predict the onset of algal blooms is clearly an issue of 
great practical importance. 

Several mechanisms of both biological and physical 
nature contribute to blooms, not all of them probably 
known yet. From the point of view of biology, the domi- 
nant controlling factor for phytoplankton growth is usu- 
ally considered to be grazing by zooplankton [14] (though 
nutrient transport, light, and the presence of a thermo- 
cline surely play an important role [15]). This is the so 
called mismatch issue [16]: since the life cycle of phyto- 
plankton is typically one order of magnitude shorter than 
that of zooplankton, a positive fluctuation in phytoplank- 
ton productivity is not compensated by a simultaneous 
increase in the zooplankton population. This produces 
the bloom event: the phytoplankton population escapes 



zooplankton control and grows to the carrying capacity 
of the medium, before also the zooplankton population 
grows appreciably, and is able to bring that of the phy- 
toplankton back to its pre-bloom level. 

At a sufficiently coarse grained scale, the system dy- 
namics can be described by phytoplankton and zooplank- 
ton concentration fields P(x,t) and Z(x.,t) [17, 18]. This 
constitutes a mean field approximation for the stochas- 
tic individual dynamics; at different levels of complexity, 
such a description may include the nutrient and even the 
detritus concentration fields N(x,t) and D(x,t) (NPZ 
[19, 20] and NPZD [21] models). The tininess of both 
phyto- and zooplankton individuals, as well as the size of 
typical areas of interest for the study of the concentration 
dynamics (at least several meters), suggests that a mean 
field description is indeed the most appropriate. The fact 
that bloom events may be triggered by fluctuations in the 
system, however, should suggest caution. 

The point that we want to examine in the present pa- 
per is precisely whether microscopic fluctuations, disre- 
garded in a mean field approach, can resurface at macro- 
scopic scale and contribute to trigger bloom events. This 
is an example of demographic fluctuation induced break- 
up of mean field descriptions of reaction-diffusion sys- 
tems, a type of phenomenon that have received a great 
deal of attention in recent years [22-26] 

In the present paper we shall consider a simple PZ 
model due to Truscott and Brindley (TB model), which 
has the nice characteristics that the plankton behaves 
like an excitable medium [18]. The original model is zero- 
dimensional, but can be easily generalized to include spa- 
tial effects such as advection and diffusion [27]. In the 
present analysis a spatially homogeneous domain will be 
considered with no advection, but with diffusive terms 
accounting for small scale motions in the water column. 
The game to play will be to determine the demographic 
fluctuations implicitly neglected in the TB model, and 
see if they can destabilize the system even at the level of 
spatial averages. 

This paper is organized as follows. In Sec. II, the 
main properties of the TB model are reviewed. In Sec. 
Ill, the expression for the stochastic contribution to de- 
mography are derived. In Sections IV and V, the effect 
of demographic noise on the dynamics of the TB model, 
without and with seasonal forcing, is analyzed. Section 
VI is devoted to the conclusions. Some technical details 



2 



on the master equation treatment of demographic fluc- 
tuations, in spatially extended domains, are provided in 
the Appendix for reference. 



II. THE TRUSCOTT-BRINDLEY MODEL 

We review here the main properties of the TB model. 
For the moment we disregard the spatial structure of the 
fields and focus on the original zero-dimensional version 
of the model, which is described by the equations: 



P = r a p(l-^)-R n 
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The main characteristics of the dynamics are the follow- 
ing: 

• Logistic reproductive behavior of the phytoplank- 
ton, with the carrying capacity K determining the 
maximum concentration that the medium can sup- 
port at steady state. 

• Grazing by zooplankton characterized by a so 
called Holling-III kind of behavior [28]. Zooplank- 
ton are able to graze on phytoplankton with opti- 
mal rate R m Z, only if the concentration of the sec- 
ond is above the level fixed by the half-saturation 
concentration a. Below this threshold, the grazing 
rate is quadratic in P: R m (P/a) 2 Z, reflecting both 
a reduced grazing ability of the zooplankton in a di- 
lute environment, and the actual reduced amount 
of food available. 

• A carrying capacity of the medium supposed much 
larger than the half saturation concentration, K 

a, meaning that at high values of P, the ability 
of zooplankton to control phytoplankton growth is 
limited. 

• A zooplankton reproductive dynamics supposed 
slower than that of the phytoplankton: /J./R m "C 1, 
while r /i? rn ~ 1. A conversion efficiency 7 as- 
sumed consistently small. 

The default values of the constants that are utilized in 
the zero-dimensional case arc [18]: 



ro = 0.3/day, R rl 
K = 108mgC/m 3 , 
7 = 0.05, 



a 



0.7/day, fi = 0.012/day, 
5.7mgC/m 3 , (2) 



where the units "mg C" stand for carbon milligrams in 
dry weight. From here we can extract three independent 
dimcnsionless groups 



Rm 

a 



0.43, 



jR-n 



0.34, 



e = — ~ 0.053. 



(3) 




FIG. 1: Trajectories approaching the fixed point (Pf,Zf) 
(fat dot to lower left corner of picture) starting from differ- 
ent initial conditions; going from outer to inner trajectory: 
(P,Z) = (0.5,0.5); (0.5.0.6); (0.6,0.66); (0.6,0.68). The pa- 
rameters are those of Eq. (2) . The dotted line indicate change 
of signature of the Jacobian matrix for Eq. (1), and could 
roughly be identified as the point where the bloom begins. 



We see that the dynamics is characterized by two inde- 
pendent small parameters: 7 and e. The dynamics de- 
scribed by Eq. (1) has a fixed point at the scale of the 
half-saturation constant a: 
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For small 7,6 C 1 and q < 1/2, this fixed point is glob- 
ally attracting (a Holling-III functional form for grazing 
appears to be crucial for stability). However, if the ini- 
tial zooplankton concentration is too low, before reaching 
the fixed point, the system will make an excursion to the 
high P ~ K range, which could be interpreted as a bloom 
event. The situation is illustrated in Fig. 1. As shown 
in figure, the onset of bloom could roughly be identified 
in the PZ plane by the line where the largest eigenvalue 
of the Jacobian of Eq. (1) crosses to positive, and phase 
points start to separate exponentially. 

The above picture of bloom triggering by zooplankton 
depletion can be improved including the effect of seasonal 
forcing. Model equations (1) can accommodate this ef- 
fect by letting the phytoplankton productivity r become 
dependent on the temperature, as suggested in [29]. The 
parameterization that we adopt is the same as in [29] : a 
Van't Hoff kind of dependence for r [30] : 

r -> r(T) = r 2^ T - T °> (5) 

and a sinusoidal dependence on time of the temperature: 

T{t) =T + ATsm(Qt + (f>). (6) 

with Q = 27r/(365 days) to allow for an annual cycle, and 
4>/{2-k) = 0.59, to have that setting t — on January 
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FIG. 2: Bloom and no-bloom cycles (heavy and thin lines 
respectively) in the seasonally forced TB model. The phase 
points circle counterclockwise along the cycle. Notice that the 
intersections of the two cycles corresponds to phese points 
that in the two cycles are associated to different (although 
close) instants of time. The dotted vertical line gives the posi- 
tion of the fixed points that would correspond to the different 
values taken by f during the year. The fat dot still identifies 
the fixed point corresponding to r = ro- The parameters in 
the graph are those in Eqs. (5-6). 



FIG. 3: Dependence on the critical forcing amplitude AT cr i t 
on the parameters 7 (thin line), e (heavy line) and q (insert). 
In the three cases AT 7jejq indicate the value of AT cri t ob- 
tained varying 7, e, q respectively and keeping the remaining 
parameters fixed to the values in Eq. (3). 



ical phytoplanktcr. A reasonable estimate (with large 
variations) for the mass of a copepod could be, for in- 
stance: 



1st, causes the first temperature minimum to occur on 
March 1st and the first maximum on August 29th. For 
the sake of defmitenes, as in we set AT = 6°C and 
v = 0.1 °CT 1 . 

Adding a seasonal forcing, turns out to modify the dy- 
namics in important way, with the single fixed point in 
the autonomous case leaving way to two stable limit cy- 
cles [29]. The situation is illustrated in Fig. 2: a small 
amplitude no-bloom cycle coexists with a large ampli- 
tude bloom cycle, each one characterized by a well de- 
fined (time-dependent) basin of attraction. As illustrated 
in Fig. 3, the small cycle will remain stable only if 
the seasonal temperature excursion AT is below a crit- 
ical threshold AT cr jt(ro, q, e, 7), that, for the values of 
the parameters quoted in Eqs. (2) is AT crit ~ 6.1°C. 
Similarly, it can be shown that, if AT is too small, only 
the no-bloom cycle will survive. The no-bloom destabi- 
lization threshold is very close to the value of the forcing 
AT = 6°C considered in [29]. The analysis in that paper 
showed in fact that addition of a fluctuating component 
to the forcing, lead to random switch from year to year 
between the two regimes. As it is clear from Fig. 2, the 
destabilization is likely to take place near the intersection 
between the bloom and no-bloom trajectories, where the 
separation between the phase points of the two cycles (at 
equal times) is smaller. 

III. DEMOGRAPHIC FLUCTUATIONS 

We consider a situation in which the typical size of a 
zooplankton individual is much larger than that of a typ- 



m z ~20 M gC, (7) 

corresponding to a typical individual size in the millime- 
ter range [6, 31]. In a no-bloom regime P,Z ~ a with a 
given as in Eq. (2), Eq. (7) would lead to a numerical 
density of the order of one zooplankton individual per 
liter. Given the much smaller size of microscopic algae, 
phytoplankton could in turn be treated as a continuum 
at those scales. We thus expect that the demographic 
fluctuations in the system be driven by the zooplankton. 

Locally, demographic fluctuations will be the result of a 
competition between stochasticity at the individual level 
of the birth-death process, and mixing by spatial trans- 
port. We shall consider a two-dimensional situation, in 
which the vertical structure of the water column is not 
resolved, and parameterize horizontal mixing through a 
diffusivity n, supposed equal for both phyto- and zoo- 
plankton. This is probably the only viable strategy to 
describe a very complex situation, in which microscopic 
swimming, stirring by larger organisms, turbulence gen- 
erated by perturbations at the water surface, all play an 
important role [32]. Likewise, we shall neglect all effects 
of large scale advection, including the forcing induced by 
the formation of fronts, where the two plankton popula- 
tions may get out of balance [5] . 

Including the possibility of a 2D spatial structure, the 
original model equations (1) can be written in the form 

P = B P (P,Z)- D p (P,Z) + kV 2 P, 

2 = B z (P,Z)-D z (P,Z) + k,V 2 Z, (8) 

where Bpz and Dpz give the local birth and death rates 
in the population, and P = P(x,t) and Z = Z(x,t) give 
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the plankton concentration averaged over the height of 
the water column. Expressing time in units i?" 1 and 
concentrations in units a (and lengths in terms of some 
reference scale), the parameters entering Eq. (8) can be 
written in the form 

B P = rP; D P = erP 2 +P 2 Z/(l + P 2 ); 

B z = 7 P 2 Z/(1 + P 2 ); D z = iq Z. (9) 

From now on, unless otherwise stated, all relations will 
be expressed in dimensionless form. 

Rather than working in a field theoretical setting [33] , 
we prefer to apply the standard system size expansion 
approach of van Kampen directly to the master equation 
for the system [34]. Let us therefore partition the do- 
main in volumes of horizontal size Ax, each one contain- 
ing instantaneously Np,Nz individuals of the P and Z 
groups. We can introduce instantaneous concentrations 
coarse grained at horizontal scale Ax, 

with 

n P , z = h(Ax) 2 /mp Z (10) 

parameterizing the size of the population in the volumes, 
and h indicating the height of the water column. Indicate 
with P = P — P and Z = Z — Z the fluctuating part of 
the coarse-grained field; in the present situation of fluc- 
tuations driven by zooplankton discreteness, we expect 
P/P,Z/Z = 0(N Z 1/2 ). Let us define normalized fluc- 
tuation fields (j), C with this scaling contribution factored 
out: 

n p = n P p :=n P (p + n z 1/2 (/)), 
n z = n z z :=n z (z + o z 1/2 c). 

The standard approach is to hypothesize that the birth 
and death rates Bp^z and Dp y z at the population scale, 
reflect birth and death rates at the individual level, con- 
ditioned to the instantaneous values of the concentration 
fields P(x.,t) and Z(x,t). At the population level, the 
transition probabilities dV in an interval dt will be: 

N P -> N P + 1 : dV = n P B P (P, Z)dt, 

Np — > Np — 1 : dV = fl P Dp(P, Z)dt, 

Nz^Nz + 1: dV = n z B z (P, Z)dt, 

Nz^Nz-l : dP = n z D z (P,Z)dt. (11) 

A standard procedure [34] leads then, from Eq. (11), to 
the master equation for the probability density function 
(PDF) p({4>i, d}, t) (the index i labels the volumes in the 



domain): 




where fi c = £l z fip, and the additional term ILp 
accounts for the exchange of plankton between adja- 
cent volumes produced by diffusion (sec Appendix). Ex- 
panding to lowest order in Op X z , both the exponen- 
tials and the reaction rates in Eq. (12) [recall that 
Bp = B P (P + n z 1/2 <f>, Z + n z 1/2 (), . . .], we obtain the 
Fokker-Planck equation: 




where ^ = [B z (Pi, Zi)+D z (Pi, Zi)]6ij, and the a PP , . . . 
give the entries of the Jacobian matrix for Eq. (1) 

(a P p = dP/dP, ...). Notice that the only noise corre- 
lator present in the equation, S^ , is the one associated 
with the variable £, i.e. with the zooplankton fluctua- 
tions. 

Going back to the original variables P, Z, and tak- 
ing the continuous limit, we see that the Fokker-Planck 
equation (13) is equivalent to the system of Langevin 
equations 

dP 

— + a PP P + a PZ Z = kV 2 P 
at 

dZ 

— + a Z pP + a zz Z = K V 2 Z + t (14) 

where 

(£(x,i)£(0,0)>=H(x,t)<S(x)J(t), 
E(x,t)=m z [B z (P,Z)+D z {P,Z)], (15) 

and we have put rhz = m z /h, P = P(x,t) and Z = 
Z(x,t) (more details in the Appendix). 

Equations (13-15) provide a mesoscopic description for 
the system dynamics that is valid as long as the fluctu- 
ations at scale Ax can be considered small. It is imme- 
diately clear that this condition cannot be satisfied in a 
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bloom regime. In fact, as illustrated in Fig. 1, bloom 
is associated with change of signature in the Jacobian 
matrix of Eq. (1) and with exponential growth of fluctu- 
ations in Eq. (14). We shall deal with this regime in Sec. 
V; let us concentrate for the moment the case of small 
fluctuations. 



IV. DYNAMICS NEAR THE FIXED POINT 

Let us consider first the case of a TB model without 
seasonal forcing, and study the fluctuations around the 
fixed point given by Eq. (4). The analysis is similar 
to the one carried on in [26] on another PZ model (the 
Levin-Segel model [35]), which focused on the destabi- 
lization of Turing patterns. The evolution equations for 
the correlation functions Cpp(x, t) — (P(x, t)P(0, t)), . . . 
can be obtained from Eqs. (14-15), and take the form in 
Fourier space: 

TjCpPk + ( fl pp + K^ 2 )Cppk + a PZ Cpzk = 

Cpzk + a ZP Cppk + (a PP + a zz + 2nk 2 )Cpzk 

+a PZ Czzw = 
^Czzk + cizpCpzk + (a zz + nk 2 )Czzk = 5. (16) 

At the fixed point we find immediately, using Eq. (9) 
and working to lowest order in e and 7: 



At small scales, the fluctuations are smeared out by dif- 



CLpp 
Q Z P 



-f(l-2q), a PZ = -q, 
2^(1 -q), a zz = 0, 

S ~ 2jm Z y/q/(l-q). (17) 

Equation (16) becomes at steady state: 

(a PP + nk 2 )Cppk + cLp Z Cpzk = 0, 

a Z pCppk + (a PP + 2nk 2 )Cpzk + CLp Z Czzk = 0, 

a ZP Cpzv. + Kk 2 Czzk = S, 

that has solution 



Cpp k ~ A-^S, 

Cpzy ^ A~ 1 a PZ (-a PP + nk 2 )E, 

Czz* A-^-app + nk 2 ) 2 -, (18) 



with 



A = appa,p Z a Z p + a pp K,k —2a PP Kk +nk. (19) 

From Eqs. (18-19), we see that there is a long wavelength 
range, dominated by demography; using Eq. (17): 



C 



ppk ^ 



qm z 



C 



f 2 (l-g)(l-2g)Vl-g' 
rh z 



PZk 



Czzw ^ 



fo(i - q) V 1 - 

rh z {l - 2 <7) / T 



g(l-g) \jl-q 



(20) 



fusion, with the asymptotic law Cppk ~ appE/ (nk 2 ) 3 , 
Cpzk — a,p Z E/(nk 2 ) 2 , Cpzk — S/(kA; 2 ). The transition 
occurs at nk 2 ~ a Z pa PZ /ap P , which sets the crossover 
length, from Eqs. (3) and (17), back to dimensional units: 

A c = V^Tm- (21) 

This is the typical distance travelled by a zooplankter in a 
lifetime and corresponds to the characteristic wavelength 
of the chemical waves supported by the system in the 
mean field. Notice that, contrary to the case considered 
in [27], for the choice of parameters utilized, no Turing 
instability is present. 

The correlation spectrum that has been obtained, char- 
acterized by a plateau at k\ c < 1, and a decay at k\ c > 1, 
corresponds to fluctuations with a correlation scale A c . 
The fluctuation amplitude can be estimated approximat- 
ing the decay at k\ c > 1 with a step function and ap- 
proximating the solution for k\ c < 1 with Eq. (20). This 
gives for the fluctuation amplitude 



Cpp(0)= J^ijCppk] ~ A- 2 C 



PPO, 



with Cppo as given in Eq. (20), and similar expressions 
for Cp Z (0) and C Z z{0). From Eqs. (20) we get for the 
ratio of the fluctuation amplitude to the mean (back to 
dimensional units): 



Cpp(O) C PZ (0) C ZZ (0) m z 



P 2 



PZ 



Z 2 



ahX 2 



(22) 



that is the ratio of the zooplankter mass and the typical 
total plankton mass in a water column of height h and 
horizontal extension A c . 



V. DESTABILIZATION OF THE NO-BLOOM 
REGIME 

Let us pass to consider a seasonally forced situation 
and ask whether demographic noise is able to destabilize 
the global bloom and no-bloom cycles. 

As already discussed, the linearized theory of Sec. IV 
cannot be utilized in the present the trajecto- 

ries evolve for most of their time in an unstable region. 
Thus, numerical solution of the master equation (12) 
becomes necessary, either by direct means, or by Mon- 
tecarlo techniques. In alternative, one may try to ex- 
tend the Langevin equation (14) beyond the linear regime 
by replacing the Jacobian matrix with the full RHS 
(right hand side) of the mean field TB model (8). It 
is convenient to express lengths in units A c , so that the 
forced equation can be written in the form 



P = Bp{P,Z) 
Z = B Z {P,Z) 



Dp(P,Z) 
D Z (P,Z) 



- AgV 2 P, 

XqV 2 Z + C (23) 



[Notice the dependence of the terms on RHS of Eq. (23) 
on the fluctuating fields P and Z; similarly, the depen- 
dence on P and Z in Eq. (15) is replaced by one on P 
and Z\. 
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FIG. 4: No-bloom cycle destabilization threshold in func- 
tion of AT, from numerical simulation in a periodic domain 
200AX x 200AX with AX = A c /10. The no-bloom cycle 
is considered destabilized if the system crosses the threshold 
P = 10 before t = 10 years. All parameters except AT and 
rhz set to the values in Eqs. (2) and (7). Initial conditions 
set equal to (Pf, Zf) uniformly in the domain. 



In this way, the role of effective diffusivity is played by 
the product Xq, while the parameter rhz, which, through 
Eq. (15), determines the amplitude of the noise £, takes 
the form, in terms of dimensional parameters: 



mz 



rrizfJ. 
ahn 



(24) 



This quantity will play the role of a control parameter 
for the theory. [We notice by the way that rhz coincides 
with the amplitude ratio in Eq. (22)]. 

The approximation leading to Eq. (23) is clearly un- 
controlled, as no account is taken of the modifications 
produced in the noise term. Nevertheless, except for very 
large values of rhz, Montecarlo simulation with Eq. (11), 
and direct solution of Eq. (23), have given very close re- 
sults. This appears to be due the fact that the noise 
plays a crucial role in destabilizing the cycles, only in 
the narrow regions around the "intersections" (see Fig. 
2). In the unstable regions, the dynamics is dominated 
by amplification of the separation between phase points 
produced by the deterministic part of the equation; the 
noise plays in the unstable regions a minor role. Thus, 
the fact that in these regions the noise in Eq. (23) is not 
well approximated, does not produce dangerous effects. 

As expected, a sufficiently high level of noise destabi- 
lizes the small cycle and leads to locking the system in the 
large bloom cycle. As illustrated in Fig. 4, the thresh- 
old in rhz becomes lower as the critical forcing AT cr i t is 
approached. 

For default values of the parameters, such as those in 
Eqs. (2) and (7), with AT = 6°C, the threshold would be 
at rhz — 4 • 10~ 5 , which, for a depth h ~ 5 m, would cor- 
respond to a diffusivity n ~ 0.2m 2 /day = 0.024cm 2 /s 
and a correlation length A c ~ 4.1m. In comparison, 
K ~ 0.01cm 2 /s would be the diffusivity that would be 
produced by microscopic swimming at speed ~ lmm/s 
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FIG. 5: Snapshots of the phytoplankton concentration field 
from a simulation with rhz = 10~ 4 ; AT = 6°C and other 
parameters set as in Eqs. (2) and (7). Size of the do- 
main 200AX x 200AX; periodic boundary conditions, with 
AX = 0.3A C . Initial conditions set equal to (Pf,Zf) uni- 
formly in the domain. A: day 700 (December 1st: no-bloom 
condition); B: day 890 (May 9th: pre-bloom condition); C: 
day 924 (June 13th: bloom peak); D: day 970 (July 29: con- 
centration minimum after bloom). 
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FIG. 6: Same as Fig. 5 in the case of the zooplankton field. 



with a persistence time between change of directions of 
the order of one second [6]. 

As illustrated in Fig. 5, even when the system is locked 
globally on a bloom cycle, its dynamics is characterized 
by spatial fluctuations, with regions of size ~ A c in which, 
during bloom events, P remains well below its typical 
bloom value. 

This picture is confirmed by looking at the evolution 
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FIG. 7: Top figure: evolution of the mean concentra- 
tion fields P (thin line) and Z (heavy line). Bottom fig- 
ure: evolution of the normalized RMS fluctuations op — 
(P 2 ) 1/2 /P (thin line), a z = {Z 2 ) 1/2 /Z (heavy line), c PZ = 
sign({PZ})\{PZ)/{PZ)\ 1/2 (dotted line). Same choice of pa- 
rameters as in Figs. 5 and 6. 




destabilization could be just a finite size effect, and would 
disappear in the thermodynamic limit. One suggestion 
that this is not the case comes from the invasive char- 
acter of the destabilization phenomenon. The sequence 
B — C in Fig. 5 gives a hint of the process: regions of size 
~ A c , characterized by high values of P, through diffusive 
coupling, destabilize nearby regions with low values of P, 
and push them in the bloom cycle. This picture is cor- 
roborated by the behavior of the system in the noiseless 
rhz = case. As shown in Fig. 8, in the absence of noise, 
the phenomenon persists: phytoplankton transfer across 
a ~ A c distance, from a bloom to a no-bloom region, 
destabilizes the second one. Thus, if a bloom bubble is 
sufficiently large (on the scale of A c ), it will survive diffu- 
sive transfer and gradually expand at the expenses of the 
surroundings. It is possible to see that the situation is 
confirmed in the case of a system with just two homoge- 
neous compartmens coupled diffusively, one evolving on 
a bloom cycle, the other on a no-bloom cycle: the no- 
bloom cycle will always be destabilized and the system 
will lock on a global bloom cycle. 
This suggests the following picture: 

• Demographic noise continuously generates local 
fluctuations in P and Z (especially in Z) that may 
push some regions of space in a bloom regime. 

• If the noise amplitude is sufficiently high, some of 
these regions will be large enough for not being 
destroyed at once by diffusivity (see passage from 
year 1 to year 2 in Fig. 8). 

• At this point diffusion, through phytoplankton 
transfer, quickly destabilizes the surrounding no- 
bloom regions. Noise is not expected to be impor- 
tant any more in this phase. 



FIG. 8: Evolution due to the effect of diffusion, in the ab- 
sence of noise, of domains characterized by bloom dynamics 
(yellow in figure; a bloom event at a given pixel is identified 
by crossing during the year of the threshold P — 12). The 
initial condition at t — was a random distribution in space 
of values P, Z in the bloom and no-bloom basin of attraction. 
Initial fraction of bloom points: 0.235 (for larger fractions, 
the system locks immediately on a bloom-cycle; the opposite 
for lower fractions). Domain characteristics and parameters 
as in Figs. 5-6. 



of the RMS fluctuations, as depicted in Fig. 7. The 
stronger fluctuation level in the phytoplankton concen- 
tration field, at the onset of the bloom events and soon 
after their disappearance, that can be seen in case B and 
D in Fig. 5, is paralleled by the double peaks in dp in 
Fig. 7. 

The above picture of global destabilization of the no- 
bloom cycle is based on numerical evidence from sim- 
ulations in a finite domain. One may question whether 



VI. CONCLUSIONS 

One can imagine a hierarchy in the description of a 
plankton population, that goes, increasing the degree 
in refinement, from zero-dimensional models, to models 
based on the dynamics of concentration fields, such as 
P and Z 1 to individual based models (IBM), taking full 
account of demographic stochasticity. In some ways, PZ 
models (and alike) may be seen as the mean field approx- 
imation of some IBM. The correspondence between the 
two levels of description is established imposing that the 
individual reproduction probabilities in the IBM (condi- 
tional to the instantaneous local values of the concen- 
tration fields), have the same form as the corresponding 
mean rates at the population level. 

By definition, a mean field approximation breaks-up 
when fluctuations are important. Usually, this means 
large amplitude fluctuations. If the mean-field dynamics 
is close to an instability, however, this is not necessarily 
the case, and global changes in the system can be induced 
by small fluctuations. 
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We have provided an example of such a situation in the 
case of a PZ model (the TB model [18]), with destabi- 
lization by demographic noise resulting in switching the 
system globally from a no-bloom to a bloom seasonal cy- 
cle. [The parameter determining the strength of the noise 
is the the dimcnsionless constant rhz of Eq. (24)]. From 
the point of view of the PZ model, this corresponds to a 
decrease of the instability threshold in the seasonal tem- 
perature forcing with respect to the mean-field case (see 
Fig. 4). For shallow water conditions (few meters depth), 
and mixing in the water column produced by diffusion 
(diffusivity in the range 0.1m 2 /day), demographic noise 
appears sufficient to cause switching between regimes. 
The effect is of the same magnitude as that of the global 
temperature fluctuations considered in [29] (fluctuation 
amplitude ~ 1 °C with correlation time equal to 35 days). 

As in other systems [25, 26], the modifications of the 
mean-field dynamics, produced by demographic noise, 
are expected to be maintained in the case of an infi- 
nite systems (thermodynamic limit). The reason is partly 
trivial: the reproduction rates at the population level are 
nonlinear functions of the concentration fields P and Z, 
and fluctuations lead necessarily to their renormalization. 
Such a picture, however, is incomplete, as the destabiliza- 
tion process, unless the demographic noise is very large, 
destabilizes the system only locally. [The characteristic 
scale of the fluctuations coincides with that of the Tur- 
ing patterns of the system, and has nothing to do with 
stochastic demography]. Only later, by an invasive pro- 
cess, which appears to be insensitive to system size, the 
destabilized regions inglobate the inactive surroundings, 
leading to a global bloom cycle. 
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Appendix A: Master equation treatment of diffusion 

We consider for simplicity a one-dimensional domain 
and a single species, say Z. Discretize the domain and 
indicate with JVj the number of individuals in slot i. We 
can define the PDF's: 

pn(n) ee PN ({nz i + n 1 ' 2 t i ,i = i,...K}) 

Suppose that individuals are transferred diffusively from 
slot i to slot i ± 1 with a rate Wi-^+i = Wi-n-i = 
Wi = Niki, where hi = Ki/(Ax) 2 , with «j = k(x^) the 
diffusivity and Ax the width of the slot. The master 



equation for pn can be written in the form 

£_ ?{ £«p{«rv. (5 £.-£) } 

x W t+k -2W l } PN . (Al) 
To derive an equation for p = p^ , we exploit the relation 
n K/2dm = dp_ n i /2 y2, d P 



Substituting into Eq. (Al) and expanding to 0(0, 1 ), 
we obtain: 



dt 4^ dQ 



d_ 



1/9 d \ 2 i I 

+ { 



d _l 

x (^"I^K 1 ^}}^ (A2) 

where we have introduced the rate density Wi = 
fi _1 Wj = Ziki. Equation (A2) could be further sim- 
plified exploiting the relation 

^[w l+ i + - 2w l ]p = 0. 

i 

At this point we write u>i = u>i + Vt^^hiQ, Wi = kiZi, 
and expand Eq. (A2) in powers of Q,. We find, to 

oin- 1 / 2 )-. 

[2i - {w i+1 + Wi-! - 2wi)}^ = 0, (A3) 

which gives, taking the continuous limit Ax — > 0, the 
diffusion equation 



dZ(x,t) _ d 2 (K(x)Z(x,t)) 



(A4) 



dt dx 2 
Going to next order, we find the equation for the fluctu- 
ations 



dp 

dt 



(A5) 



where 



11, = Ki-i— Q-i + Ki+i— Q+i 

oQ-i oQi+i 
i~ \ ® 

- (Ki-lCi-l + K i+lQ+l)-^7- 
+ Ki+lA+ll— 



+ 



d(i 



(A6) 
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From Eqs. (A5-A6) we obtain the equation for the cor- 
relation C y (t) = <C*(*)Cj(*)>: 

Cjk = — 2(kj + hk)Cjk + kj-iCj-i^k + Kk-iCj^k-i 



+ Kj+lCj+l t k + Kt+lCji+l + -j^jk 



(A7) 



where 



7Tjfc = < 



{wj-i +2wj +w j+ i}, k=j, 
-{%_!+%}, k=j-l, 
0, fc < j - 1, 



and 7Tjfe = 7rfej . Taking the continuous limit we find that 
the source term Wij is infinitesimal: 

TTij -> -4AxZ(x)k(x)(5"(2; - j/) ->• 0, 



and the remaining terms give the diffusion equation for 



dC(x,y;t) d 2 ( K (x)C(x,y;t)) 



dt 



+ 



dx 2 

d 2 ( K (y)C(x,y;t)) 
dy 2 



(A8) 



The corresponding dynamics for the fluctuation field £ is 
diffusive as well: 



d(( X ,t) d 2 (K(x)C(x,t)) 



dt 



dx 2 



(A9) 



which leads to the diffusion terms to the RHS of the 
Langevin equations (14). 



[1] C.B. Field, M.J. Beherenfeld, J.T. Randerson, and P.G. 

Falkowski, Science 281, 237 (1998) [21 
P.J.S. Franks, Limnol. Oceanogr. 42, 1297 (1997) [22 
N. Blackburn, T. Fenchel, and J. Mitchell, Science 282, 
2254 (1998) [23 
A.P. Martin, Philos. Trans. R. Soc. London A 365, 2873 
(2005) [24 
R. Reigada, R.M. Hillary, M.A. Bees, J.M. Sancho and 

F. Sagues, Proc. R. Soc. Lond. B 270, 875 (2003) [25 
A.M. Metcalfe, T.J. Pedley and T.F. Thingstad, J. Ma- 
rine Sys. 49, 105 (2004) [26 
M. Levy, Lect. Notes Phys. 774, 219 (2008) 
A. Bracco, S. Clayton and C. Pasquero, J. Geophys. Res. [27 
114, C02001 (2009) 

W.J. McKiver and Z. Neufeld, Phys. Rev. E 83, 016303 [28 
(2011) 

P.J.S. Franks, Limnol. Oceanogr. 42, 1273 (1997) [29 
U. Siegenhalter and J.L. Sarmiento, Nature 356, 589 

(1993) [30 
T.J. Smayda, Limnol. Oceanogr. 42, 1137 (1997) 
OS. Yentsch, B.E. Lapointe, N. Poulton and D.A. Phin- [31 
ney, Harmful Algae, 7, 817 (2008) 

M. Scheffer, S. Rinaldi, Y.A. Kutznetsov and E.H. Van [32 
Nes, Oikos 80, 519 (1997) 

T.J. Smayda, Limnol. Oceanogr. 42, 1132 (1997) 
D.H. Cushing, Adv. Mar. Biol. 26, 249 (1990) [33 

G. T. Evans and J.S. Parslow, Biol. Oceanogr. 3, 327 
(1985) 

J.E. Truscott and J. Brindley, Bull. Math. Biol. 56, 981 

(1994) [34 
J.H. Steele and E.W. Henderson, Am. Naturalist 117, 
676 (1981) [35 
M. Fasham, H. Ducklow and S. McKelvie, J. Marine Res. 



48, 591 (1990) 

A.M. Edwards, J. Plankton Res. 23, 386 (2001) 

WR. Young, A.J. Roberts and G. Stuhne, Nature 412, 

328 (2001). 

A.J. McKane and T.J. Newman, Phys. Rev. Lett. 94, 
218102 (2005) 

CR. Doering, K.V. Sargsyan and L.M. Sander, Multi- 
scale Modelling and Simulation 3, 283 (2005) 
T. Butler and D. Reynolds, Phys. Rev. E 79, 032901 
(2009) 

T. Butler and N. Goldenfeld, Phys. Rev. E 84, 011112 
(2011) 

I. Siekmann and H. Malchow, Math. Model. Nat. Phe- 
nom. 3, 114 (2008) 

C.S. Holling, Can. Entomol. 91, 293 (1959); C.S. Holling, 
Can. Entomol. 91, 385 (1959) 

J.A. Freund, S. Mieruch, B. Scholze, K. Wiltshire and U. 
Feudel, Ecol. Complex. 3, 129 (2006) 
J.A. Berges, D.E. Varela and P.J. Harrison, Mar. Ecol. 
Prog. Ser. 225, 139 (2002) 

E.J. Gonzalez, T. Matsumura-Tundisi and J.G. Tundisi, 
Braz. J. Biol. 68, 69 (2008) 

To be precise, in order for the two species to have iden- 
tical diffusion properties, we could assume that passive 
transport be dominant over swimming. 
M. Doi, J. Phys. A 9, 1465 (1976); A.S. Mikhailov, Phys. 
Lett. 85, 214 (1981); N. Goldenfeld, J. Phys. A 17, 2807 
(1985); L. Peliti, PJ. Physique 46, 1469 (1985); H.K 
Janssen and U.C. Tauber, Ann. Phys. 315, 147 (2005) 
N.G. Van Kampen, Stochastic processes in physics and 
chemistry (Elsevier, New York, 1992) 
S.A. Levin and L.A. Segel, Nature 259, 659 (1976) 



